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We study stability conditions of the full Hamiltonian constraint equation describing the quantum 
dynamics of the diagonal Bianchi I model in the context of LQC. Our analysis has shown robust 
evidence of an instability in the explicit implementation of the difference equation, implying impor- 
tant consequences for the correspondence between the full LQG theory and LQC. As a result, one 
may question the choice of the quantisation approach, the model of lattice refinement, and/or the 
role of the ambiguity parameters; all these should in principle be dictated by the full LQG theory. 



^ ■ I. INTRODUCTION 

° : 

Loop Quantum Gravity (LQG) [1[ is a non-perturbative, background independent, canonical quantisation of General 
Relativity in four space-time dimensions. Even though the full theory of LQG is not yet complete, its successes 
encourage the application of LQG techniques to mini-superspaces obtained by a symmetry reduction. The application 
I ^ , of LQG to the cosmological sector is known as Loop Quantum Cosmology (LQC) @, in the homogeneous and 
isotropic cosmological models with a massless scalar field, which plays the role of an internal time parameter according 
which other physical quantities "evolve" , quantum geometry effects of the full LQG theory lead to a repulsive force 
in the Planckian regime. Thus, the big bang singularity is resolved and replaced by a quantum bounce [4]. The 
r— 1| underlying discreteness of LQC is the key element for the existence of the quantum bounce; similar results have thus 
^ been also obtained in the context of other models. 

I ' LQC quantum dynamics are determined by a difference, rather than a differential, equation, as a result of quantum 
£h , geometry effects. However, such effects can be neglected as one departs from the Planckian regime, and quantum 
i dynamics can then be well approximated by the Wheeler-DeWitt (WDW) differential equation. LQC is formulated 

in terms of SU(2) holonomics of the connection and triads. In the "old" quantisation, the quantised holonomies were 
taken to be shift operators with a fixed magnitude, but later it was found that this leads to problematic instabilities in 
the continuum semi-classical limit, where the WDW wave-function becomes a good approximation to the difference 
equation of LQC. In a dynamical equation closer to what is expected to be obtained from the full LQG theory, lattice 
refinement would take place during the evolution, since full Hamiltonian constraint operators generally create new 
vertices of a lattice state in addition to changing their edge labels. The effect of the refinement of the discrete lattice 
has been modelled and the elimination of the instabilities in the continuum era has been explicitly shown [1, 0] . 
Lattice refinement leads to new dynamical difference equations which, in general, do not have a uniform step-size 
making their study quite involved. In contrast to isotropic models, which can be understood in terms of wave-functions 
on a one-dimensional discrete mini-superspace, anisotropic models with higher-dimensional mini-superspaces, can be 
• • ■ more subtle. For the partial difference equations of anisotropic models, stability issues can turn out to be more serious 
than in isotropic ones, leading to consistency tests, and thus restricting possible quantisation freedom. In Ref. Q we 
have proposed a numerical method, based on Taylor expansions, which provides the necessary information to calculate 
the wave-function at any given lattice point. We have developed Q numerical schemes for both the one-dimensional 
homogeneous and isotropic cosmological case, which has analytic solutions, as well as the two-dimensional case of a 
Schwarszchild interior, which cannot be exactly solved. 

LQC issues of the Bianchi type I models, the simplest among anisotropic cosmologies, have been also investigated. 
Besides their simplicity, such models are very interesting for addressing the issue of space-like singularities in the 
context of the full LQG theory. As in the isotropic case, a massless scalar field plays the role of an internal time 
parameter. Recent analysis @ has shown that the big bang singularity is solved by quantum gravity effects, while 
LQC dynamics is well approximated by that of the WDW theory once quantum geometry effects become negligible. 

The aim of this paper is to analyse the stability conditions of the solutions to the full Hamiltonian constraint. 
Unstable (i.e., growing) solutions would indicate unphysical spurious solutions, for which there is no correspondence 
between the difference (valid in the LQC regime) and the differential (WDW) equations. This would indicate an 
inconsistency between the full LQG theory and the mini-superspace LQC approach, implying the possibility of a 
weakness of the employed quantisation approach. This work is organised as follows: In Section II we outline the basic 
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formalism of LQG and LQC. In Section III we perform a stability analysis. We summarise our results and we discuss 
the outcome of our findings in Section IV. 



II. BASICS OF THE LQG/LQC FORMALISM 

Let us restrict ourselves to diagonal Bianchi I metrics, for which space-time metric in Cartesian coordinates, r, X4 
(i=l,2,3), reads 

3 

ds 2 = -N 2 dr 2 + a f dx l > (!) 

i=l 

where N is the lapse function and a; (with i = 1, 2, 3) stand for the three directional scale factors. Following Ref. @ 
we choose r to satisfy Or = 0. 

LQG/LQC are based on a Hamiltonian formulation of General Relativity, with basic variables an SU(2) valued 
connection A l a and the conjugate momentum variable which is a densitised triad E? , a derivative operator quantised 
in the full LQG theory in the form of fluxes. As for any quantisation scheme based on a Hamiltonian framework 
or an action principle, for the homogeneous flat model one should regularise the divergences which appear due to 
the homogeneity as the action and Hamiltonian are integrated over spatial hyper-surfaces. We thus restrict spatial 
homogeneity and Hamiltonian to an elementary cell V, which we choose so that its edges lie along the fixed coordinate 
axis Xi (with i = 1, 2, 3). In addition, we fix a fiducial flat metric °q a b, with line element 

3 

ds 2 =J2^ 2 i- (2) 

i=i 

The lengths of the three edges of the elementary cell V and its volume, as measured by the fiducial flat metric °q a b, 
are denoted by Li (with i — 1,2,3) and Vo = L1L2L3, respectively. 

The densitised triad carries information about the spatial geometry, encoded in the three-metric, while the connec- 
tion carries information about the spatial curvature, in the form of the spin-connection and the extrinsic curvature. 
We introduce physical triads °ef and their dual (°e™ a uj 3 a — Sj) co-triads °uj l a — D a x l , satisfying °q a b = °uj\ 
Note that i refers to the Lie algebra index and a is a spatial index with a, i = 1, 2, 3. The physical co-triads are given 
by uf — a 1 °to l a , and the physical three-metric by q a b = 

The six-dimensional phase space is defined through the SU(2) connection A l a and the triad Ef given by 

E? = piLiV- 1 ^^ , (3) 

where the connection components Ci and the momenta pi are constants; q — (P1P2P3) qV^ 1 stands for the determinant 
of the physical spatial metric q a b- The three momenta pi are related to the three scale factors through 

pi = sgn(ai)|a 2 a 3 |L 2 i3 
p 2 = sgn(a 2 )|aia 3 |LiL 3 

p 3 = sgn(a 3 )|aia2|iii 2 • (4) 
The pairs c l ,pi (with i = 1,2, 3) satisfy the Poisson brackets relations: 

{ C >,} = 87rG 7 <5} , (5) 

with 7 the Barbero-Immirzi parameter. 

Two of the constraints of the full LQG theory, namely the Gauss and the diffeomorphism constraints are identically 
satisfied and one is therefore left with the Hamiltonian constraint, as for the isotropic case. Restricting the integration 
to the fiducial cell V, the Hamiltonian constraint reads 

C = I N(n grav + H mat tcr)d 3 x , (6) 
Jv 

where 7i gra v and 7i m attcr stand for the gravitational and the matter parts of the constraint densities, respectively. 
The lapse function N is iV = \f\piP2~P~s\- 
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Since Bianchi I models are spatially flat, the matter part of the Hamiltonian constraint can be written as @ 



Hgrav = ~ „ 2 r T , {PlP2ClC 2 +P1P3C1C3 +P2P3C2C3) . (7) 

8nGj z ^pip2P3Vo 
The matter part of the Hamiltonian constraint is Q 

matter 1 (8) 

where /3 ma ttor is the matter energy density of the matter field, chosen to be a massless scalar field T; 

2 

Pt 

Pmattcr — ^1 f > (9) 

with the canonically conjugate momentum of T. The scalar field T can be considered as an evolution parameter 
in the classical theory, and as a viable internal time parameter in the subsequent quantum theory. The justification 
for this choice lies in the fact that since px is a constant of motion, T grows linearly in time r, for any solution to the 
field equations. The full Hamiltonian constraint, Eq. ([7]) can then be finally written as Q 



"H = ~ 8ttG , 7 2 ( piP2ClC2 +P1P3C1C3 +P2P3C2C3) + ^ . (10) 



Let us proceed with the quantum kinematics of Bianchi I LQC. The gravitational part of the kinematic Hilbcrt 
space, Wf[„ v , can be expressed in the momentum, pi (with i = 1,2,3), representation. Given an orthonormal basis 
states \pi,P2,P3), which are eigenstates of quantum geometry, consider a linear combination 



with finite norm, namely 



and 



I*) = *(Pl>P2,P3)bl,P2,P3> , (11) 



\V(pi,P2,P3)\ 2 <oo , (12) 



(13) 



The action of the elementary operators, which are the three momenta pi (with i — 1,2,3) and the holonomies 
along edges parallel to the three axis a;, (with i = 1,2,3) — completely determined by almost periodic (£ is any real 
number) functions exp(i£cj) of the connection — is given by @ 

Pl\Pl,P2,P3) = Pl\Pl,P2,P3) 

exp(Ud)\pi,p2,P3) = \pi - 8-kGjM,P2,P3) , (14) 

and similarly for p2,exp(z£c2) and P3, expiiic^). 

One has then to build the quantum analogue of the Hamiltonian constraint, along the lines of the isotropic case. 
To do so, one has to find the operator on the gravitational sector of the kinematic Hilbert space, corresponding to 
the curvature F ab k of the connection A l a , given by 

F ab k = 2d [a A b] k + e tJ k AlAl . (15) 

As it is known from the isotropic case, the connection operator does not exist in LQG/LQC; we cannot take the 
limit of the area enclosed by a plaquette to go to zero, since the minimum area enclosed by the plaquette is the 
nonzero eigenvalue App l (with A a dimensionless number, A = 4-^/3^7) of the area operator. To single out a unique 
plaquette of the many ones enclosing an area Appj on each of the three faces of the elementary cell V, we will use 
the natural gauge fixing available for the diagonal Bianchi I case, and a correspondence between kinematic states in 
LQG and LQC. In this way, one obtains that the curvature operator reads Q 

p k_ f kf sin M c 0, , V f sin M c 0, , V nfn 
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where 



with 



sm/ic o \ sm//c l 



r u a = — ^— "ul , (17) 



Mi = 



M2 = 



'NAZ 2 



pi 



|j*2P3 I 



'| P2 |AP 



pi 



\PiPd,\ 



The functional dependence of /Etj on is essential since otherwise quantum dynamics can depend on the choice of the 
fiducial cell V. 

Consequently, one can now write the quantum analogue of the full Hamiltonian constraint, Eq. It reads 0] 

- h 2 d^{\, T) = 9#(A, T) , (19) 

where = — C grav . 

To simplify the gravitational sector of the Hamiltonian constraint, one can introduce the volume of the elementary 
cell V as one of the arguments of the wave function. Let us then set Q 

v = 2AiA 2 A 3 , (20) 

which is directly related to the volume of V, namely 



F*(Ai, A a , v) = 2tt| 7 | VA>|^*(Ai, A 2 , v) , (21) 

with 7 = sgn(piP2P3)|7|- Thus, the new configuration variables will be Ai, A 2 , v. 

In the next section, we will write out explicitly the full Hamiltonian constraint and we will then study the stability 
of its solutions. 

III. STABILITY ANALYSIS 

The basic difference equation arising from the loop quantisation of the Bianchi I model reads Q 

vrG 



$*(Ai,A 2 ,i/;T) = — v^[(v + 2)v^+4*i"(Ai,A2,i/;T)-(i/ + 2)v^*3"(A 1 ,A 2 ,i/;T) 



- (i/ - 2) >/J7*o (Ai,A a ,i/;T) + (i/-2)vlI^*l(Ai,A2,i/;T)J , (22) 

where 

*+(Ai,A 2 ,^;T) = J2 ^Mi.ijAj,^^) 

i#i=(0,l,2) 

*4 (A 1 ,A 2 ,i/;T) = ^ *(o i A 1 ,a J A 2 ,i/-4;T) 

i#i=(-3,-2,0) 

*+(Ai,A 2 ,^;T) = VfaXuajX^viT) 

^=(-1,0,1) 

*o (Ai,A 2 ,i/;T) = ^ *(a,:Ai,a,A 2 ,^;T) , (23) 

#i=(-2,0,3) 
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and the functions a* have been defined as follows: 



v — 4\ ( v — 1\ ( v 



a-3 = ^ , a-2 = , a-i 



v-2) ' \ v ) ' V^ + 2 

ao-1, «^(^)' ° 2 -(^)' fl3 =(^)' (24) 

Numerical evolution can in principle be carried out by restricting to the positive octant (Ai > 0, A 2 > 0, v > 0), thus 
eliminating the sgn(Aj) factors which are otherwise appearing in various terms. 

Here we wish to examine the stability of the vacuum solutions, in which case the solution is static, namely 
* ( Ai , A 2 , v\ T) = f ( X 1 , A 2 , v) , and Eq. becomes 



* 4 + (Ax, A,, u) = (Ax, A,, u) + (^f ) (A x , A 2 , v) - (^f) ") , (25) 

for v ^ 0; otherwise the above equation must be multiplied by y^, thus corresponding to the classical singularity. 
The geometry of this difference equation is drawn in Fig. [TJ Equation ([25]) can be used to evaluate the value of the 
wave-function on the u + A plane, given suitable boundary conditions on the v and v — A planes. The requirement that 
the arguments must be positive (i.e., Ai > 0, A 2 > 0, v > 0) reduces the required number of boundary conditions. For 
the purpose of our work, it is sufficient to consider starting from a plane in which v — A > 0. 

In addition to specifying the boundary conditions on the v and v — A planes, we are also required to specify the value 
at five of the points given in f + (Ai,A a ,i/). There are in total 23 values that are required and with such initial data 
the difference equation, Eq. I|22p . can be used to evaluate the 24 th point. Once this point has been evaluated, it can 
be used to "move" the central point and evaluate the wave-function at subsequent positions in the v + 4 plane. In this 
way the difference equation can be used to find the wave- function that is consistent with the Hamiltonian constraint, 
Eq. (f2"2"l) . and the boundary conditions. In principle, this procedure can be iterated to evaluate the consistent wave- 
function for all subsequent i/-planes, however the stability of the difference equation can be investigated even at this 
first iteration. 

As shown in Fig. [TJ there is a choice to be made as to which point in the v + 4 plane is to be calculated from the 
difference equation. This choice amounts to deciding whether to increase Ai or A 2 first, when populating the v + 4 
plane. From the point of view of the plane, the difference equation, Eq. I|22p. can be seen as progressively evaluating 
the wave-function at points first along either the Ai direction or the A 2 one (see, Fig. [2]). In this sense, we can consider 
Eq. (|22[) as an "evolution" equation of a wave- function with respect to either Ai or A 2 , subject to suitable boundary 
conditions. It is important to realise however that this "evolution" has only to do with the order in which the points 
are evaluated and is not related, in any way, to evolution of the wave-function with respect to time. 

With this view, standard von Neumann stability analysis can be preformed on Eq. (|22|) . to see if the system is 
stable [HHHI- Here however caution is necessary. Von Neumann's analysis is typically used to see if there are growing 
mode solutions to a particular discretised version of an underlying differential equation. In this case, the difference 
equation is the fundamental evolution equation, which can be approximated by a differential equation (the anisotropic 
Wheeler-DeWitt equation Q) in a suitable limit. In standard numerical implementations of differential equations, 
the stability of the system is important only because artificial numerical rounding errors can grow to dominate the 
behaviour of the solution, however the situation here is very different. In principle, the difference equation, Eq. (I22p . 
is exact and hence all solutions should be considered, however in practise we wish to restrict only to those solutions 
that closely approximate General Relativity at large scales. This makes the use of von Neumann stability analysis 
useful, since we are comparing a particular difference equation, with the differential equation it approximates, however 
it is important to remember that the motivation is very different than in standard numerical analysis. 

For homogeneous and isotropic cosmologies, a local stability analysis of the corresponding difference equation to 
determine the behaviour of spurious solutions was performed in Ref. [l(J, using higher order spin J representations 
of the holonomies for the quantisation. It was found [T(| that the use of higher spin holonomies to regulate the 
gravitational part of the constraint operator leads to modifications, which are qualitatively similar to those of the 
inverse scale factor. Stability analysis has shown that the J = 1 difference equation is not locally stable. To further 
determine whether these spurious solutions represent a problem with the quantisation, the authors of Ref. [l(| have 
studied the physical inner product, since unphysical solutions would have either vanishing or infinite physical norm 
and would be modded out of the physical Hilbert space. For the cases of Bianchi I locally rotationally symmetric 
cosmology and that of the Schwarzschild interior geometry, a von Neumann stability analysis of a difference equation 
obtained by a previous quantisation approach was carried out in Ref. [e|, where there were identified large regions 
in space-time that have generically instabilities. In what follows, we will look for spurious solutions to Eq. (|^)) . in 
the sense that they do not approximate solutions to the relevant Wheeler-DeWitt equation in the large volume limit. 
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FIG. 1: The geometry of the points used in the difference equation that results from the Hamiltonian constraint, for the 
Bianchi I model. 



As in standard von Neumann stability analysis, we will decompose the solutions of the difference equation, Eq. 125)1 . 
into Fourier modes and look for growing modes. Specifically, we consider the ansatz 

* (Ai, A 2 , v) = T (Ar) exp (t (u;X 2 + X i/)) , (26) 

where we have chosen the Ai direction to be the direction in which the v + 4 plane is "evolved" . Using the above 
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FIG. 2: The difference equation gives us a point in the v + 4 plane, given the required 23 points. Exactly which point is 
calculated via the difference equation is somewhat arbitrary and essentially describes the way in which the v + 4 plane is 
calculated. In the l.h.s. scheme (a) the point (02 Ai, 0,1X2, v + 4) is calculated, in which case the v + 4 plane is evaluated first 
along constant A2. In the r.h.s. scheme (b) the point chosen is (aiAi, 02 A2, v + 4) and the v + 4 plane would be evaluated first 
along constant Ai. 



ansatz, Eq. (|25|) becomes 



ijij=(0,X,2) ' V+ i#j=(-l,0,l) 



7 v i#j=(-2,0,3) 



7 .W=(-3,-2,0) 

(27) 

To simplify each of the summations, we proceed as follows: 

T( ai Ai)e lwA2 = T(aiAi) (e^ 2 + e iWB - lA ") 

i^j=(-l,0,l) 



+ T(a Ai) (e 

+ T (a_iAi) ( e <wolA2 + e luA2 ) , (28) 



-ia;a_iA2 j_ gicjaiA2 
ia;aiA2 r A 



which becomes 



^ T(a,i\\)< 

^'=(-1,0,1) 



iujX'2 



2e 



T(aiAi)e 5 C os I '- 

+T (a Ai) cos (uj\ 2 [a x - 1)) 

icA 2 (a 1 -i) ) /wA 2 (ai-l) 
+T(aiAi)e 2 cos' 



where we made use that 

a_ 3 - 1 = - (a 3 - 1) , a_ 2 - 1 = - (a 2 - 1) , o_i - 1 = - (oi - 1) . 

We can simplify the other summations in a similar way. 

Explicitly putting in the values of oi, a 2 , 03 given in Eq. (|24p . the difference equation, Eq. (|27p . becomes 



T a A! e cos 7\T' )+T{a 1 X 1 ) e~ cos — 



">-1 I uj\ 2 

+T(a 2 Ai)e-+2 cos 

v + 2 



1/ + 4 



z/-2 



T (a_!Ax) e^ 2 cos f + T (oqAx) cos ^ 2 ^ 2 



+T(aiAi) e "+ 2 cos 



z/ + 2 

ljA 2 
i/ + 2 



T (a_ 2 A x ) e^f cos ( + T (a Ar) e^fy cos f ^ 



- 2 



■"->>2 / (jjA 2 

+i(a3Ai)e " cos 

V v 



v{v-2) 



v-2\ \v-A\ 



v + 2 J V z/ + 4 



- 'j / wA 2 \ „ , , , -i»*2 / wA 2 

T a_ 3 Ai cos — - + T (a_ 2 Ai) e^ - cos - 

z/ / \ z/ — 2 



+T(a Ai)e -<- 2 ) cos' 



V i^(^-2) 

Up to this point the equation is exact, however expanding in terms of small 1/z/, Eq. (|3 1 [) becomes 

e(4 X ) 



2 

1 

z/ 



1 - - \e {-A X ) 
v 



T (o Ai) e (2A) cos (2A) + (T (01A2) + T (a 2 Ai)) e (A) cos (A) 
T (a Ai) cos (2A) + T (a_iAi) e (-A) cos (-A) + T (01A1) e (A) cos (A) 
1 - T (a_ 2 Ai) e (-A) cos (-A) + T (a 3 Ai) e (A) cos (A) 
T (a_ 3 Ai) e (-A) cos (-A) + T (a_ 2 Ai) e (-A) cos (-A) + T (a Ai) e (-A) cos (-A) 

-o (v ' 



where we have defined the function 



e (x) = e 1 ' 



and the variable 



A = LUX2/V . 

Equation (j32|) can be re-ordered to read 

AT (a 3 Xx) = BT (a 2 A x ) + CT (oiA x ) + DT (a Ai) + ET (o_iAi) + i^T (a_ 2 Ai) + GT (o_ 3 Ai) 

where 



.4 
B 
C 

D 

E 
F 
G 



I- 6 - 

v 



: (A) cos (A) 



-e (4x) e (A) cos (A) 



1 - - - e (4 X ) 



(A) cos (A) 



-e (4 X ) e (2A) + 1 - - 
v 



cos (2A) 



1 

1/ 



2" 
1 

v 

I- 6 - 

V 



e (-A) cos (-A) 



e (-4%) 



e (-A) cos (-A) 



e(-4x)e(-A) cos(-A) 



Equation 



is equivalent to the vector equation 



where we have defined the vectors 



AhT 3 - M 2 T 2 , 



T(ajAi) 
T (aj-iAi) 
T(aj_ 2 Ai) 
T(aj_ 3 Ai) 
T(ai_ 4 Ai) 
r(ai_ 5 A!) 



for i = 2, 3 



and the matrices 



M, 



/A 













1 

















1 

















1 

















1 





V 











l) 



M 9 



I B C D E F G\ 

1 

1 

1 

1 

V 1 J 

Stability of this system is then given by the eigenvalues of the matrix (Mi) -1 M 2 . More particularly, if 

max I A I < 1 V lo and x 1 

where A are the eigenvalues of the matrix (Mi)" 1 M 2 , then the amplitude T (c^Ai) is less than that of previous 
namely the difference equation is stable. 
One finds, in block form, that 
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where I5 is the 5x5 identity matrix, 5 is the zero vector and 



A = 

B = 

C = 

D = - 

E = - 



1 



e(4 X ) 



e(4x) 



1 



e (Ax) e (2A) + 1 - 



e (-2A) 



6 



e(-A) 



1 <l-->e (-A X ) 



cos (2A) 
cos (A) 



e (-2A) 



1 - 



e(-4 X ) 



F = 



1 



e(-4 X )e(-2A) , 



(42) 



with e(x) as defined in Eq. (f3"3")> . previously. The eigenvalues of Eq. (|4"T]) are found by solving the characteristic 
equation 



(Mi) 'Ms-Ale 



= 0, 



(43) 



for the eigenvalues A; note that 1q is the 6x6 identity matrix. We are looking for the maximum |A|, for all u and 
X- We can immediately see that the system will not be stable, since the inverse of Mi only exists when \A\ ^ 0. The 
cases when \A\ = correspond to 



A = 



(2n- 1)tt 



or, equivalently, using Eq. 



u> 



(2n - 1)tt v 
2 AT 



(44) 



(45) 



with HgZ and these modes are explicitly unstable. This can be understood by noting that the amplitude T(<Z3Ai) is 
multiplied by A, which can be made arbitrarily small, hence then the amplitude T (a 3 Ai) has to be arbitrarily large. 

We can go further and consider the th order limit in the (l/^) — > expansion, in which the definitions given in 
Eq. g2D simplify to 



fj(0 



— g4ix 



= (e itx+2lA + 1) e 



;2A 



cos A 



-4ix 



-2iA 



— e -4ix-2'iA 



(46) 



where the superscript (°\ reminds us that we are working to the th order in the small (1/^) expansion. 

If we further consider the modes given by A = 7r/4 and x = 0, then the above coefficients, Eq. (I46p . become simply 



i(o) = i i 5(0) = o 
= 1 , P = i 

= , = - 



(47) 
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In this specific case, the matrix given in Eq. (|4Tj) reads 



M- 1 M 2 



( 1 1 i -i\ 

1 

1 

1 

1 

Vo 1 / 



(48) 



the determinant of which is simply 
implying 



det (Mf x M 2 



(49) 



(50) 



Xj are the eigenvalues of the matrix M x x Mi. Since |n| =1 A,-| = 1, either max f|Aj|J > 1, or |Aj| = 1, Mj with j = 

1, • • • ,6. We can rule out the second possibility by explicit evaluation of the characteristic equation, Eq. (|43|) . for this 
ansatz. 

To be more specific, set Xj = exp(idj) and solve Eq. (|43)l . subject to the limit (1/v) — > 0, for the modes A = 7r/4 
and x = 0, to find 9j. In this case, Eq. (|43| becomes 



i9j \ „§iQj 



i (1 - e 2 '^') =0 , Vj with j = 1, • • -6 



(51) 



The above equation, Eq. (|5ip . has only two (numeric) solutions, which without loss of generality we denote by 6\, 82, 
and are approximately equal to 6\ = 1.18123 and 82 = 2.30716, for 6\, 82 £ (0, 2ir). However, using Eq. (|50[) . the sum 
of the phases of the six eigenvalues must satisfy 



^2 °j = (2n - 1) 7r for neZ 

3=1 



(52) 



With only two solutions, the eigenvalues must be degenerate. Let us suppose that |Ay| = 1, Vj with j = 1, • • • ,6, and 
consider p eigenvalues with phase #1 and q eigenvalues with phase #2, where p and g are integers satisfying p + q = 6. 
We can then look for any combination of degeneracies (i.e., any values of p and q) that satisfy Eq. (|52")l . Explicitly 
it can be verified that there is no such solution, which implies that not all of the eigenvalues lie on the complex unit 
circle and hence there must be at least one eigenvalue with |Aj| > 1. 

A partial pro of of this result in the general case can be produced by using a variant of the Gershgorin circle 
theorem [l3l. Il4|. The standard theorem states that the eigenvalues of a matrix M. = ((%), lie within the i discs, 
D (an, R) (called Gershgorin discs) in the complex plane with centre an and radius R = Ylj-fj \ a ij\- It cari further be 
shown that if the discs are disjoint, then there is at least one eigenvalue within each connected region. For the case 
of the matrix given by Eq. (|4"Tj) this implies that all of the eigenvalues lie within the discs 



^(0,1) 



D 



(a,\ 



B\ 



I CI 



\D\ + \E\ 



\F 



(53) 



Of the two discs, the second one is the most interesting. It is centred at A and one can easily check that for 7^ 0, 

it is beyond the unit complex circle, i.e., \A\ > 1. However, one can also check that the radius satisfies 



\B\ + \C\ + \D\ + \E\ + \F\ > \A\ - 1 , 



(54) 



except for small values of v. Thus, the two Gershgorin discs intersect and we cannot say that there is an eigenvalue 
with |Aj| > 1. However, by noting that \C\ becomes arbitrarily large for A — > tt/2, one realises that the radius of 
the second disc in Eq. (|53|) . encompasses all of the complex plane. This would tend to suggest that there is at least 
one eigenvalue that is not constrained to have |Aj| < 1. A variation on the proof of the standard Gershgorin circle 
theorem can be used to show that this is indeed the case. 

Consider the case of a matrix M. = (a^) such that |ai3 1 ~> J2j^3 \ a ij\- Then the characteristic equation is 



E 

3=1 



(55) 
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where x = (xi) is the eigenvector of M. and A is the corresponding eigenvalue. Expanding this sum as 



JV3 



(56) 



gives 



A - a.,3 — 



^3 



which is valid, provided Xi ^ 0. If we take Xi to be 

Xi = max(xj) for j ^ 3 

we have 



A - a i3 



■E3 



< E K 

j¥3 



(57) 



(58) 



(59) 



Thus, the eigenvalue A is within a disc, centred at the point aax^/xi with radius given by the sum of the magnitudes 
of the elements along the i th row of Ai, excluding the third element. In particular, if jaia] f oo, then for x^/xi > 0, 
the centre of the disc tends to infinity. Provided the sum Ylj^ remains finite, the eigenvalue A will lie within 
a disc that is entirely outside the complex unit circle and hence |A| > 1. This is precisely the situation we have for 
M = M^ 1 M 2 , in the case of A — > 7r/2. 

The final element that is required for this proof is that x^/xi > or, more precisely, that 013X3/0:1 3> Y] ,y 3 |oij|, 
given that | ai3 1 3> Ylj-f^ \ a ij \ ■ ln the particular case of the matrix given by Eq. (141[) . we can evaluate the simultaneous 
equations implied by the characteristic equation, Eq. (|43p . to find 



C'x 3 
X3 



Xxi 

A24 



Xi 
Xi 



Xx 2 



X2 

x 5 



Xx 3 
Xx 6 



where we have used the approximation that C dominates the terms in J^- |oii| 

1/3 



This gives 



^3 

aid, — 

Xi 



C 



(60) 



(61) 



Thus, provided that 



C 



1/3 



»E 



3^3 l Ul il 



the proof is valid and we have max(|Aj|) > 1. Note that this condition is 



certainly met as A — > tt/2, since C diverges, whilst |oii| remains finite. This is essentially the result we preempted 
in the comments following Eq. (|4"5)) . however here we have explicitly extended it to the case of C large, but not infinite 
(i.e., the case when Mi is invertible, but A is large). 



IV. CONCLUSIONS 



The aim of this paper is to study the stability of the Hamiltonian constraint equation valid for anisotropic Bianchi I 
LQC. Performing a von Neumann stability analysis, we have shown that if the difference equation admits solutions 
with amplitudes that grow locally, then it is not locally stable. On the one hand, this result certainly questions the 
validity of the quantisation, since any semi-classical solutions would quickly become dominated by the expanding 
spurious ones. On the other hand however, the presence of such an instability may not be, necessarily, a problem, 
since it might be that the unstable trajectories are explicitly removed by the physical inner product. 

More precisely, the difference equation, given by Eq. (f2"2"]) , is unconditionally unstable. By this we mean that there 
is no region of (Ai,A2,^) in which the difference equation, Eq. (|22[) . is stable. It is worth noting however, that in 
Eq. (|55|) we choose to re-order the difference equation in such a way that it produces a single amplitude (T(o3Ai) 
in Eq. (|35p ). given the other 23 amplitudes. This is clearly an explicit implementation of the equation. It is also 
possible that this difference equation could be implemented via an implicit scheme, i.e., that the equation could be 
re-ordered to give (say) two amplitudes, given the values of the other 23 or 22 amplitudes. In order for the system to 
give solutions, one would then have to implement consistency relations between the calculated amplitudes at different 
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iterations. There are, of course, many ways that such an implicit implementation of the difference equation could be 
under taken and they could, in principle, have different stability properties. 

We have demonstrated the presence of an instability in the explicit implementation of the difference equation, 
Eq. (|22j) . in several ways: we have first shown that for a particular set of critical modes, A = (2n — 1)tt/2 1 the system 
is unstable. We have then showed that in the large v limit, the system is again unstable for the modes A = 7r/4 and 
X = 0. Finally, we have formally showed that the system is unstable for a general v, for modes that approach the 
critical value. This was done via a version of the Gershgorin circle theorem, which have explicitly demonstrated the 
instability, even for modes approaching (but not reaching) the critical value. 
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